package scu.maqiang.possion;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T41;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.*;

public class Possion3DT41 {

	public Possion3DT41(Mesh3T4 mesh) {
		this.mesh = mesh;
		fs = new FES3T41(mesh);
		A = new SRMatrix(fs.GetNdof());
		x = new double[fs.GetNdof()];
		gradient = new double[3][fs.GetNdof()];
		RHS = new double[fs.GetNdof()];
	}
	
	
	public static void main(String[] args) {
		
		/*
		Mesh3T4Generator mg = new Mesh3T4Generator();
		Mesh3T4 mesh3D = mg.cube(50, 50, 50);
		GeneralCoefFunc ConstFunc = (xyz, relabel, param) -> 1.0e0;
        GeneralCoefFunc sourceFunc = (xyz, relabel, param) ->  
                          2 * xyz[0] * (1 - xyz[0]) * xyz[1] * (1 - xyz[1]) 
                        + 2 * xyz[1] * (1 - xyz[1]) * xyz[2] * (1 - xyz[2])
                        + 2 * xyz[2] * (1 - xyz[2]) * xyz[0] * (1 - xyz[0]);
		Possion3DT41 p2DT31 = new Possion3DT41(mesh3D);
		p2DT31.DirichletProblem(BVPType.COMMON, ConstFunc, null, sourceFunc, null, null, null, null, "MyResult3D.dat");
		*/

		/*
		int N = 8;
		Mesh3T4Generator  mgt3 = new Mesh3T4Generator();
		Mesh3T4 fineMesh = mgt3.cube(40, 40, 40);
		fineMesh.saveMesh("Cube.msh");
		ScalarFunc regionFineFunc = (xy, llabel, param) -> {
            double xCell = xy[0] * N - (int)(xy[0] * N);
            double yCell = xy[1] * N - (int)(xy[1] * N);
            double zCell = xy[2] * N - (int)(xy[2] * N);
            return (Math.abs(xCell - 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25 && Math.abs(zCell - 0.5) < 0.25)? 1.0:-1.0;
        };
        fineMesh.setDomainLabel(regionFineFunc, null, 10);
        
        ScalarFunc cf = (xyz, label, param) -> label == 10? 0.01: 1;
        ScalarFunc sourceFunc = (xy, llabel, param) -> 10.0;
        Possion3DT41 p2dt31 = new Possion3DT41(fineMesh);
        p2dt31.DirichletProblem(BVPType.COMMON, cf, null, sourceFunc, null,
        		                null, null, null, "MyNewFineSolution3D.dat");

		 */

//		Mesh3T4Generator mgt4 = new Mesh3T4Generator();
//		Mesh3T4 cubeMesh = mgt4.readMesh("Cube.msh", MeshFormat.msh);
//		cubeMesh.toTecplot("MyCube.dat");
		Mesh3T4 mesh3D = new Mesh3T4().cube(16, 16, 16);
		ScalarFunc ConstFunc = (xyz, relabel, param) -> 1.0e0;
//		GeneralCoefFunc sourceFunc = (xyz, relabel, param) ->
//				2 * xyz[0] * (1 - xyz[0]) * xyz[1] * (1 - xyz[1])
//						+ 2 * xyz[1] * (1 - xyz[1]) * xyz[2] * (1 - xyz[2])
//						+ 2 * xyz[2] * (1 - xyz[2]) * xyz[0] * (1 - xyz[0]);
		Possion3DT41 p2DT31 = new Possion3DT41(mesh3D);
		p2DT31.DirichletProblem(BVPType.COMMON, ConstFunc, null, ConstFunc, null, null, null, null, "MyResult3D.dat");
//		Mesh3T4 holeMesh = mgt4.readMesh("HoleCompo512.msh", MeshFormat.msh);
//		holeMesh.saveMesh_mesh("HoleMesh.mesh");

//        FES3T41 fs = new FES3T41(mesh3D);
//        SRMatrix A = new SRMatrix(fs.GetNdof());
//
//        GeneralCoefFunc regionCoef = (xyz, relabel, param) -> relabel == label ? 0.01 : 1.0;
//        GeneralCoefFunc ConstCoef = (xyz, relabel, param) -> 1.0e0;
//        GeneralCoefFunc[] gcf = {regionCoef};
//        gcf[0] = ConstCoef;
//        
//        GeneralCoefFunc coefFunc = (xyz, relabel, param) -> 
//                    2 * xyz[0] * (1 - xyz[0]) * xyz[1] * (1 - xyz[1])
//                  + 2 * xyz[1] * (1 - xyz[1]) * xyz[2] * (1 - xyz[2])
//                  + 2 * xyz[2] * (1 - xyz[2]) * xyz[0] * (1 - xyz[0]);
//        fs.assembleHeatStiff(gcf, null, BVPType.COMMON, A);
//        //fs.assembleHeatStiff(new double[] {1.0}, BVPType.COMMON, A);
//        //System.out.println(A);
////        System.out.println(A.makeEntrySymmetric());
//        double[] RHS = new double[fs.GetNdof()];
//        //GeneralCoefFunc ConstCoef = (xyz, relabel, param) -> 1.0;
//        gcf[0] = coefFunc;
//        fs.assembleHeatSource(gcf, null, BVPType.COMMON, RHS);
//        A.makeEntrySymmetric();
//        fs.applyBCZero(A, RHS, Direct.X, BCType.MBN, 1, 2, 3, 4, 5, 6);
//        
////        System.out.println(A);
//
////        IterSSolver solver = new IterSSolver(1.0e-15, 1000000);
////        
//        double[] x = new double[fs.GetNdof()];
////        CSRMatrix mat = new CSRMatrix(A);
////        solver.PCGSSOR(A, RHS, x, 1.5e0, 1);
////        Arrays.fill(x, 0.0);
////        solver.PCGSSOR_2(A, RHS, x, 1.5e0, 1);
////
////        
//        //System.out.println(A);
////        
////        //System.out.println(mat);
////        
////        Arrays.fill(x, 0.0);
////        solver.CG(mat, RHS, x, 1);
////        
////        Arrays.fill(x, 0.0);
////        solver.CG(A, RHS, x, 1);
//        
//        Arrays.fill(x, 0.0);
//        NewIterSSolver newSolver = new NewIterSSolver(A, 1.0e-15, 1000000);
//        newSolver.PCGSSOR(RHS, x, 1.5e0, 1);
//
//        //Arrays.fill(x, 0.0);
//        //newSolver.GMRES(RHS, x, 200);
//        //newSolver.GMRES_SSOR(RHS, x, 1.5, 30, 2);
//        mesh3D.toTecplot("Solutions3DT4.dat", x);
//        
//        double[][] gradU = new double[3][fs.GetNdof()];
//        fs.computeGradient(x, gradU);
//        //mesh2D.toTecplot("SolutionsT32.dat", gradU);
//        System.out.println("Solution Gradient: ");
//        System.out.println("Dx min:" + MVO.min(gradU[0]) + "\t" + "max: " + MVO.max(gradU[0]));
//        System.out.println("Dy min:" + MVO.min(gradU[1]) + "\t" + "max: " + MVO.max(gradU[1]));
//        System.out.println("Dz min:" + MVO.min(gradU[2]) + "\t" + "max: " + MVO.max(gradU[2]));
//        
////        ArrayList<Integer>  li = new ArrayList<Integer>();
////        li.ensureCapacity(100);
////        for(int i = 0; i < 100; i++) {
////        	li.add(null);
////        }
////        li.add(1, 2);
////        System.out.println(li.toString());
//        
////        Mesh2DT3 mesh6 = mg.square2D(2, 2, T3Type.Left);
////        mesh6.displayBoundarys();
////        mesh6.extractBoundarysFromElements();
////        mesh6.displayBoundarys();/

	}
	
	
	/**
	 * Possion问题Dirichlet边界条件
	 * 
	 * @param kxxFunc
	 * @param paramKxx
	 * @param sourceFunc
	 * @param paramSource
	 * @param tp
	 * @param DirichletBoundaryLabel
	 * @param fixedFunc
	 * @param paramFixed
	 * @param reportFile
	 */
	public void DirichletProblem(BVPType tp, ScalarFunc kxxFunc, double[] paramKxx, ScalarFunc sourceFunc,
			double[] paramSource, int[] DirichletBoundaryLabel, ScalarFunc fixedFunc, double[] paramFixed,
			String reportFile) {
		fs.assembleStiff(new ScalarFunc[] { kxxFunc }, paramKxx, tp, A);
//		System.out.println(A.display());
//		System.exit(1);
		if (sourceFunc != null) {
			fs.assembleSource(new ScalarFunc[] { sourceFunc }, paramSource, tp, RHS);
		}

//		SRMatrix M = new SRMatrix(fs.GetNdof());
//		fs.assembleHeatMass(new double[]{1.0}, BVPType.CONSIST_MATRIX, M);
//		M.sort();
//		System.out.println(M.display());
//		System.exit(1);
		double[] MyRHS = new double[fs.GetNdof()];
		fs.assembleFlux(new double[]{1.0}, BVPType.COMMON, MyRHS, 1);
		//System.out.println(Arrays.toString(MyRHS));
		//System.exit(1);
		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);

		//System.out.println(Arrays.toString(RHS));
		if (fixedFunc == null) {
			fs.applyBC_MBN(RHS, Direct.All, 0.0, DirichletBoundaryLabel);
		} else {
			fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[] { fixedFunc }, paramFixed, DirichletBoundaryLabel);
		}
		//System.out.println("OK!");
		solveAndComputeGradientAndExport(A, RHS, x, reportFile);
	}

	/**
	 * 
	 * @param tp
	 * @param kxx
	 * @param source
	 * @param DirichletBoundaryLabel
	 * @param fixedBV
	 * @param reportFile
	 */
	public void DirichletProblem(BVPType tp, double kxx, double source, int[] DirichletBoundaryLabel,
			double fixedBV, String reportFile) {
		fs.assembleStiff(new double[] { kxx }, tp, A);
		fs.assembleSource(new double[] { source }, tp, RHS);
		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);
		fs.applyBC_MBN(RHS, Direct.X, fixedBV, DirichletBoundaryLabel);
		solveAndComputeGradientAndExport(A, RHS, x, reportFile);
	}

	public void GeneralPossionProblem(BVPType tp, ScalarFunc kxxFunc, double[] paramKxx, ScalarFunc sourceFunc,
			double[] paramSource, int[] DirichletBoundaryLabel, ScalarFunc fixedFunc, double[] paramFixed,
			int[] NeumannBoundaryLabel, ScalarFunc fluxFunc, double[] paramFlux, int[] RobinBoundaryLabel,
			ScalarFunc robinFluxFunc, double[] paramRobinFlux, ScalarFunc convectFunc,
			double[] paramConvectFunc, String reportFile) {

		fs.assembleStiff(new ScalarFunc[] { kxxFunc }, paramKxx, tp, A);
		if (sourceFunc != null) {
			fs.assembleSource(new ScalarFunc[] { sourceFunc }, paramSource, tp, RHS);
		}
		if (NeumannBoundaryLabel != null) {
			fs.assembleFlux(new ScalarFunc[] { fluxFunc }, paramFlux, tp, RHS, NeumannBoundaryLabel);
		}
		if (RobinBoundaryLabel != null) {
			fs.assembleHeatConvectSurface(new ScalarFunc[] { convectFunc }, paramConvectFunc, tp, A,
					RobinBoundaryLabel);
			fs.assembleFlux(new ScalarFunc[] { robinFluxFunc }, paramRobinFlux, tp, RHS, RobinBoundaryLabel);
		}

		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);
		fs.applyBC_MBN(RHS, Direct.X, new ScalarFunc[] { fixedFunc }, paramFixed, DirichletBoundaryLabel);
		solveAndComputeGradientAndExport(A, RHS, x, reportFile);

	}

	public void GeneralPossionProblem(BVPType tp, double kxx, double source, int[] DirichletBoundaryLabel,
			double fixedBV, int[] NeumannBoundaryLabel, double fluxValue, int[] RobinBoundaryLabel,
			double robinFluxValue, double convectValue, String reportFile) {
		fs.assembleStiff(new double[] { kxx }, tp, A);
		fs.assembleSource(new double[] { source }, tp, RHS);
		if (NeumannBoundaryLabel != null) {
			fs.assembleFlux(new double[] { fluxValue }, tp, RHS, NeumannBoundaryLabel);
		}
		if (RobinBoundaryLabel != null) {
			fs.assembleHeatConvectSurface(new double[] { convectValue }, tp, A, RobinBoundaryLabel);
			fs.assembleFlux(new double[] { robinFluxValue }, tp, RHS, RobinBoundaryLabel);
		}

		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);
		fs.applyBC_MBN(RHS, Direct.X, fixedBV, DirichletBoundaryLabel);
		solveAndComputeGradientAndExport(A, RHS, x, reportFile);
	}
	
	private void solveAndComputeGradientAndExport(SRMatrix A, double[] RHS, double[] x, String reportFile) {
		NewIterSSolver solver = new NewIterSSolver(A, true);
		solver.PCGSSSOR(RHS, x, 1.50, 1);

//		DirectSSolver solver = new DirectSSolver(A);
//		x = solver.CholeskySolve(RHS, 1);
		fs.computeGradient(x, gradient);
		System.out.println("Solution: ");
		System.out.println("x0 min:" + MVO.min(x) + "\t" + "max: " + MVO.max(x));
		System.out.println("Solution Gradient: ");
		System.out.println("Dx min:" + MVO.min(gradient[0]) + "\t" + "max: " + MVO.max(gradient[0]));
		System.out.println("Dy min:" + MVO.min(gradient[1]) + "\t" + "max: " + MVO.max(gradient[1]));
		System.out.println("Dz min:" + MVO.min(gradient[2]) + "\t" + "max: " + MVO.max(gradient[2]));
		double[][] mat = new double[][] { x, gradient[0], gradient[1], gradient[2]};
		if (reportFile != null) {
			mesh.toTecplot(reportFile, mat);
		}
	}
	
	
	public double[] exportSolution() {
		return x.clone();
	}
	
	public double[][] exportGradient() {
		double[][] gd = new double[3][];
		gd[0] = gradient[0].clone();
		gd[1] = gradient[1].clone();
		gd[2] = gradient[2].clone();
		return gd;
	}
	
	Mesh3T4 mesh;
	FES3T41 fs;
	SRMatrix A;
	double[] RHS;
	double[] x;
	double[][] gradient;
}
